home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Analog / ASPAnalyze.m next >
Encoding:
Text File  |  1992-08-18  |  13.0 KB  |  426 lines
  1. (*  :Title:    Analog Signal Processing Analysis  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    Perform basic analysis on a one-dimensional continuous-time
  7.         functions, including plotting f(t) and graphing the poles
  8.         and zeroes of F(s).
  9.  *)
  10.  
  11. (*  :Context:    SignalProcessing`Analog`ASPAnalyze`  *)
  12.  
  13. (*  :PackageVersion:  2.7  *)
  14.  
  15. (*
  16.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  17.         Georgia Tech Research Corporation
  18.  
  19.     Permission to use, copy, modify, and distribute this software
  20.     and its documentation for any purpose and without fee is
  21.     hereby granted, provided that the above copyright notice
  22.     appear in all copies and that both that copyright notice and
  23.     this permission notice appear in supporting documentation,
  24.     and that the name of the Georgia Tech Research Corporation,
  25.     Georgia Tech, or Georgia Institute of Technology not be used
  26.     in advertising or publicity pertaining to distribution of the
  27.     software without specific, written prior permission.  Georgia
  28.     Tech makes no representations about the suitability of this
  29.     software for any purpose.  It is provided "as is" without
  30.     express or implied warranty.
  31.  *)
  32.  
  33. (*  :History:    began -- March, 1990 (adapted from "ZTransform.m")    *)
  34.  
  35. (*  :Keywords:    signal processing, Laplace transform, pole-zero plot *)
  36.  
  37. (*  :Source:    *)
  38.  
  39. (*  :Warning:    *)
  40.  
  41. (*  :Mathematica Version:  1.2 or 2.0  *)
  42.  
  43. (*  :Limitation:  *)
  44.  
  45. (*  :Functions: ASPAnalyze  *)
  46.  
  47.  
  48.  
  49. (*  B E G I N     P A C K A G E  *)
  50.  
  51. BeginPackage[ "SignalProcessing`Analog`ASPAnalyze`",
  52.           "SignalProcessing`Analog`Fourier`",
  53.           "SignalProcessing`Analog`LaPlace`",
  54.           "SignalProcessing`Analog`LSupport`",
  55.           "SignalProcessing`Support`TransSupport`",
  56.           "SignalProcessing`Support`ROC`",
  57.           "SignalProcessing`Support`SigProc`",
  58.           "SignalProcessing`Support`SupCode`" ]
  59.  
  60.  
  61. If [ TrueQ[ $VersionNumber >= 2.0 ],
  62.      Off[ General::spell ];
  63.      Off[ General::spell1 ] ];
  64.  
  65.  
  66. (*  U S A G E     I N F O R M A T I O N  *)
  67.  
  68. ASPAnalyze::usage =
  69.     "ASPAnalyze[f, t] and ASPAnalyze[f, t, start, end] \
  70.     will plot f(t) as a one-dimensional, continuous-time function \
  71.     (real part shown as solid lines, imaginary part as dashed lines). \
  72.     It will also print the strip of convergence, indicate stability \
  73.     criteria, display the pole-zero diagram, and plot the magnitude \
  74.     and phase responses. \
  75.     For the two-dimensional version, t is a list of two variables \
  76.     (like {t1, t2}) and start/end are optional two-element lists \
  77.     specifying the range of the time-domain plot. \
  78.     All non-variable symbols will be assigned a value of 1 when the \
  79.     analyzer needs to generate a plot (override this by using options \
  80.     like a -> 2)."
  81.  
  82. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  83.  
  84.  
  85. Begin["`Private`"]
  86.  
  87.  
  88. (*  M E S S A G E S  *)
  89.  
  90. LaPlace::notvalid =
  91.     "The forward bilateral LaPlace transform could not be found."
  92. ASPAnalyze::notinteresting =
  93.     "Could not determine the important section of the frequency response."
  94. ASPAnalyze::hasval =
  95.     "The global variable '' has value so it could not be used \
  96.     as a transform variable; '' was used instead."
  97. CTFTransform::notvalid = "The forward Fourier transform could not be found."
  98.  
  99.  
  100. (*  ASPAnalyzeQ  *)
  101. ASPAnalyzeQ[start_, end_] :=
  102.     NumberQ[N[start]] && NumberQ[N[end]] && N[end >= start]
  103.  
  104. (*  checkvar  *)
  105. SetAttributes[checkvar, {HoldFirst}]
  106. checkvar[ xsym_, xval_, str_ ] := xsym /; SameQ[xsym, xval]
  107. checkvar[ xsym_, xval_, str_ ] :=
  108.     Block [ {sym},
  109.         sym = Unique["s"];
  110.         Message[ ASPAnalyze::hasval, str, sym ];
  111.         sym ]
  112.  
  113. (*  FourierExtent  *)
  114. FourierExtent[freqresp_, defaultvalues_, w_, mag_] :=
  115.     Block [    {newmag, wmax, wmin},
  116.         newmag = mag;
  117.         {wmin, wmax} = Extent1D[ freqresp /. defaultvalues, w ];
  118.         Which [ SameQ[wmin, -Infinity] && SameQ[wmax, Infinity],
  119.               wmin = wmax = 0,
  120.             SameQ[wmin, -Infinity],
  121.               wmin = wmax / 100,
  122.             SameQ[wmax, Infinity],
  123.               wmax = 100 wmin,
  124.             True,
  125.               newmag = Linear ];
  126.         {wmin, wmax, newmag} ]
  127.  
  128. (*  LtoFreqResp  *)
  129. LtoFreqResp[ funct_, ltrans_, defaultvalues_,
  130.          onedflag_, polezeroflag_, lrules_ ] :=
  131.     Block [ {lfreqresp, plotrange, polelist = {0}, rm, rp, stable},
  132.  
  133.         (*  Print the signal and its Laplace transform  *)
  134.         Print[ funct ];
  135.         Print[ "has the following Laplace transform:" ];
  136.         Print[ TheFunction[ltrans] ];
  137.  
  138.         (*  Print the Region of Convergence  *)
  139.         Print[ "The region of convergence is:"];
  140.         rm = GetRMinus[ltrans];
  141.         rp = GetRPlus[ltrans];
  142.         If [ onedflag,
  143.              Print[ rm, " < Re(s) < ", rp ],
  144.              Print[ rm[[1]], " < Re(s1) < ", rp[[1]] ];
  145.              Print[ rm[[2]], " < Re(s2) < ", rp[[2]] ] ];
  146.  
  147.         (*  Print stability information  *)
  148.         stable = PrintStability[ltrans];
  149.  
  150.         (*  From now on, the default values are used  *)
  151.         If [ ! EmptyQ[defaultvalues],
  152.              Print[ "The default values are now being considered." ] ];
  153.  
  154.         (*  Show the Pole-Zero Plot  *)
  155.         If [ polezeroflag,
  156.              polelist = PoleZeroPlot[ltrans /. defaultvalues] ];
  157.  
  158.         lfreqresp = TheFunction[ltrans] /. lrules;
  159.  
  160.         If [ stable /. defaultvalues,
  161.              Print[ "Since the signal is stable, the" ];
  162.              Print[ "frequency response can be computed" ];
  163.              Print[ "directly from the Laplace transform." ];
  164.              plotrange = All,
  165.  
  166.              Print[ "Since the signal is not stable, the" ];
  167.              Print[ "frequency response has no meaning, but it" ];
  168.              Print[ "will be plotted anyway.  However, the DC" ];
  169.              Print[ "(zero-frequency) term will be Infinity which" ];
  170.              Print[ "will generate a warning from Mathematica." ];
  171.              plotrange = Automatic,
  172.  
  173.              Print[ "Even though the stability of the signal" ];
  174.              Print[ "is not known, the frequency response will" ];
  175.              Print[ "still be plotted." ];
  176.              plotange = Automatic ];
  177.  
  178.         { lfreqresp, plotrange, polelist } ]
  179.  
  180. (*  formatting  *)
  181. Format[ w ]  := "w"
  182. Format[ w1 ] := "w1"
  183. Format[ w2 ] := "w2"
  184.  
  185. Format[ s ]  := "s"
  186. Format[ s1 ] := "s1"
  187. Format[ s2 ] := "s2"
  188.  
  189.  
  190. (*  PrintDefaults  *)
  191. PrintDefaults[oplist_, varlist_, defaultvalues_] :=
  192.     Block [    {},
  193.         If [ EmptyQ[oplist],
  194.              Print[ "For plotting only, these symbols will" ];
  195.              Print[ "be assigned a value of 1:  ", varlist ],
  196.              Print[ "For plotting only, these symbols will" ];
  197.              Print[ "will be assigned default values:" ];
  198.              Print[ defaultvalues ] ] ]
  199.  
  200. (*  PrintStability  *)
  201. PrintStability[ltrans_] :=
  202.     Block [    {stable},
  203.         stable = Stable[ltrans];
  204.         Which [ SameQ[stable, True],
  205.               Print[ "The system is stable." ],
  206.             SameQ[stable, False],
  207.               Print[ "The system is not stable." ],
  208.             True,
  209.               Print[ "The system is stable if ", stable ] ];
  210.         stable ]
  211.  
  212.  
  213. (*  A S P A n a l y z e  *)
  214.  
  215. ASPAnalyze[f_, t_Symbol] := ASPAnalyze[f, t, -20]
  216. ASPAnalyze[f_, t_Symbol, start_] := ASPAnalyze[f, t, start, start + 40]
  217.  
  218. ASPAnalyze[f_, t_Symbol, start_, end_, op___] :=
  219.     Block [    {bogus, complexpoles, ctft, defaultvalues = {}, freqresp,
  220.         fnumeric, ftemp, funct, ltrans, mag, oplist, plotrange = All,
  221.         polezeroflag, polelist, stable, varlist, wmax, wmin},
  222.  
  223.         funct = ToContinuous[f];
  224.         oplist = ToList[op];
  225.  
  226.         (*  CONTINUOUS-TIME DOMAIN ANALYSIS  *)
  227.  
  228.         (*  Convert expression to plottable/functional form    *)
  229.         (*  and replace constants like Pi with numerical value *)
  230.  
  231.         fnumeric = N[ TheFunction[funct] ];    
  232.  
  233.         (*  Find all valueless symbols besides the time variable *)
  234.         (*  The first parameter in Summation is a local symbol   *) 
  235.  
  236.         ftemp = fnumeric /.
  237.               ( Summation[n_, l_, u_, inc_][expr_] :>
  238.                 bogus[l, u, inc, GetVariables[expr /. n -> l]] );
  239.         varlist = Select[ GetVariables[ftemp /. oplist],
  240.                   Function[x, ! SameQ[x,t]] ];
  241.         If [ ! EmptyQ[varlist] || ! EmptyQ[oplist],
  242.              defaultvalues = oplist ~Join~ Map[(#1 -> 1)&, varlist];
  243.              PrintDefaults[ oplist, varlist, defaultvalues ];
  244.              fnumeric = fnumeric /. defaultvalues ];
  245.  
  246.         (*  Plot the signal in the time-domain  *)
  247.  
  248.         Print[ "Real part of the function is shown as solid lines.\n\
  249.             Imaginary part of the function is shown as dashed lines."];
  250.  
  251.         polezeroflag = ( ! SameQ[fnumeric, 0] );
  252.         SignalPlot[fnumeric, {t, start, end},
  253.                PlotLabel -> "Continuous-Time Domain Analysis" ];
  254.  
  255.  
  256.         (*  LAPLACE-DOMAIN ANALYSIS  *)
  257.  
  258.         (*    Find F(s); if it exists, print it out, give stability  *)
  259.         (*  information, and display a pole-zero plot.             *)
  260.  
  261.         ltrans = LaPlace[funct, t, s];
  262.         If [ ! SameQ[Head[ltrans], LTransData],
  263.              Message[ LaPlace::notvalid ];
  264.              ctft = CTFTransform[funct, t, w];
  265.              If [ ! SameQ[Head[ctft], CTFTData],
  266.               Message[CTFTransform::notvalid]; Return[ltrans] ];
  267.              freqresp = TheFunction[ctft];
  268.              wmin = wmax = 0;
  269.              mag = Linear,
  270.  
  271.              { freqresp, plotrange, polelist } =
  272.             LtoFreqResp[ funct, ltrans, defaultvalues,
  273.                      True, polezeroflag, ( s -> I w ) ];
  274.              complexpoles = Map[Im, polelist];
  275.              wmin = Min[ 0, Abs[complexpoles] ] / 10;
  276.              wmax = Max[ 0, Abs[complexpoles] ] 10;
  277.              mag = Log ];
  278.  
  279.         (*  FREQUENCY-DOMAIN ANALYSIS  *)
  280.  
  281.         (*  Plot the interesting section of the freq. response  *)
  282.         (*  That is, plot the stuff around the break points     *)
  283.  
  284.         Print[ funct ];
  285.         Print[ "has the following frequency response:" ];
  286.         Print[ freqresp ];
  287.  
  288.         If [ wmin == 0 && wmax == 0,
  289.              {wmin, wmax, mag} =
  290.                  FourierExtent[freqresp, defaultvalues, w, mag] ];
  291.         If [ wmin == 0 && wmax == 0,
  292.              Message[ ASPAnalyze::notinteresting ] ];
  293.         If [ SameQ[wmax, 0], wmax = 10 ];
  294.         If [ SameQ[wmin, 0], wmin = wmax/100 ];
  295.         MagPhasePlot[ freqresp /. defaultvalues, { w, wmin, wmax },
  296.                   Domain -> Continuous, DomainScale -> Log,
  297.                   MagRangeScale -> mag, PhaseRangeScale -> Degree,
  298.                   PlotRange -> plotrange ];
  299.  
  300.         (*  Return the Laplace transform  *)
  301.  
  302.         ltrans /. s -> checkvar[Global`s, Global`s, "s"] ] /;
  303.  
  304.     ASPAnalyzeQ[start, end]
  305.  
  306.  
  307. ASPAnalyze[f_, {t1_Symbol, t2_Symbol}] := ASPAnalyze[f, {t1, t2}, {-10, -10}]
  308. ASPAnalyze[f_, {t1_Symbol, t2_Symbol}, start_] :=
  309.     ASPAnalyze[f, {t1, t2}, start, start + {20, 20}]
  310.  
  311. ASPAnalyze[f_, {t1_Symbol, t2_Symbol}, start_, end_, op___] :=
  312.     Block [    {complexpoles, ctft, defaultvalues = {}, deltafuns, freqresp,
  313.          fnumeric, ftemp, funct, ltrans, mag, oplist, plotrange = All,
  314.          polelist = {0}, polezeroflag, slist, stable,
  315.          tlist, varlist, wmax, wmin},
  316.  
  317.         oplist = ToList[op];
  318.         tlist = { t1, t2 };
  319.         slist = { s1, s2 };
  320.         funct = ToContinuous[f];
  321.  
  322.         (*  CONTINUOUS-TIME DOMAIN ANALYSIS  *)
  323.  
  324.         (*  Convert expression to plottable/functional form    *)
  325.         (*  and replace constants like Pi with numerical value *)
  326.  
  327.         fnumeric = N[ TheFunction[funct] ];    
  328.  
  329.         (*  Find all valueless symbols besides the time variables *)
  330.         (*  The first parameter in Summation is a local symbol    *) 
  331.  
  332.         ftemp = fnumeric /.
  333.               ( Summation[n_, l_, u_, inc_][expr_] :>
  334.                 bogus[l, u, inc, GetVariables[expr /. n -> l]] );
  335.         varlist = Select[GetVariables[ftemp /. oplist],
  336.                  Function[x, ! SameQ[x,t1] && ! SameQ[x,t2]]];
  337.         If [ ! EmptyQ[varlist] || ! EmptyQ[oplist],
  338.              defaultvalues = oplist ~Join~ Map[(#1 -> 1)&, varlist];
  339.              PrintDefaults[oplist, varlist, defaultvalues];
  340.              fnumeric = fnumeric /. defaultvalues ];
  341.  
  342.         (*  Plot delta functions separately from rest of function  *)
  343.  
  344.         polezeroflag = ( ! SameQ[fnumeric, 0] );
  345.         SignalPlot[fnumeric,
  346.                {t1, start[[1]], end[[1]]},
  347.                {t2, start[[2]], end[[2]]},
  348.                PlotLabel -> "Continuous-Time Domain Analysis" ];
  349.  
  350.         (*  LAPLACE-DOMAIN ANALYSIS  *)
  351.  
  352.         (*  Find F(s) and print it out.  *)
  353.  
  354.         ltrans = LaPlace[funct, tlist, slist];
  355.         If [ ! SameQ[Head[ltrans], LTransData],
  356.              Message[ LaPlace::notvalid ];
  357.              ctft = CTFTransform[funct, tlist, {w1, w2}];
  358.              If [ ! SameQ[Head[ctft], CTFTData],
  359.               Message[CTFTransform::notvalid]; Return[ltrans] ];
  360.              freqresp = TheFunction[ctft];
  361.              mag = Linear,
  362.  
  363.              { freqresp, plotrange, polelist } =
  364.             LtoFreqResp[ funct, ltrans, defaultvalues,
  365.                      False, polezeroflag,
  366.                      { s1 -> I w1, s2 -> I w2 } ];
  367.              mag = Log ];
  368.  
  369.         (*  FREQUENCY-DOMAIN ANALYSIS  *)
  370.  
  371.         (*  Plot the interesting section of the freq. response  *)
  372.         (*  That is, plot the stuff around the break points     *)
  373.         (*  This is only possible for separable transforms.    *)
  374.  
  375.         Print[ funct ];
  376.         Print[ "has the following frequency response:" ];
  377.         Print[ freqresp ];
  378.  
  379.         complexpoles = Map[Im, Select[polelist, NumberQ]];
  380.         wmin = Min[ 0, Abs[complexpoles] ] / 10;
  381.         wmax = Max[ 0, Abs[complexpoles] ] 10;
  382.         If [ wmin == 0 && wmax == 0,
  383.              Message[ ASPAnalyze::notinteresting ] ];
  384.         If [ SameQ[wmax, 0], wmax = 10 ];
  385.         If [ SameQ[wmin, 0], wmin = 1/10 ];
  386.  
  387.         MagPhasePlot[ freqresp /. defaultvalues,
  388.                   {w1, wmin, wmax}, {w2, wmin, wmax},
  389.                   Domain -> Continuous, DomainScale -> Log,
  390.                   MagRangeScale -> mag, PhaseRangeScale -> Degree,
  391.                   PlotRange -> plotrange ];
  392.  
  393.         (*  Return the Laplace transform  *)
  394.  
  395.         ltrans /. { s1 -> checkvar[Global`s1, Global`s1, "s1"],
  396.                 s2 -> checkvar[Global`s2, Global`s2, "s2"] } ] /;
  397.  
  398.     ASPAnalyzeQ[ start[[1]], end[[1]] ] &&
  399.         ASPAnalyzeQ[ start[[2]], end[[2]] ]
  400.  
  401.  
  402. (*  A u t o c o r r e l a t i o n  *)
  403.  
  404.  
  405.  
  406. (*  E N D     P A C K A G E  *)
  407.  
  408. End[]
  409. EndPackage[]
  410.  
  411. If [ TrueQ[ $VersionNumber >= 2.0 ],
  412.      On[ General::spell ];
  413.      On[ General::spell1 ] ];
  414.  
  415.  
  416. (*  H E L P     I N F O R M A T I O N  *)
  417.  
  418. Combine[SPfunctions, { ASPAnalyze } ]
  419. Protect[ASPAnalyze]
  420.  
  421.  
  422. (*  E N D I N G     M E S S A G E  *)
  423.  
  424. Print[ "The analog signal analyzer is now loaded." ]
  425. Null
  426.